# 
# data <- readRDS('data_sim_p_1000.rds')
# for(x in names(data)) assign(x, data[[x]], envir = globalenv())
# U_full <- data$U
# rm(data)
# 
# S.data.time <- S.data$obs
# S.data.time.log.sum <- sum(log(S.data.time))
# Exemple d'optimalité entre les 3 fct MH, MH_high_dim et MH_high_dim_para sur le modèle non lineaire mixte
m <- function(t, phi1, phi2, phi3) (phi1  )/(1+exp((phi2-t)/phi3))
#=======================================#
p <- 1000
parameter <- list(sigma2 = .05^2,
                  mu = c(0.9,90,5),
                  omega2 = c(0.005, 40, 1),
                  bara = 90,
                  barb = 30,
                  baralpha = 0.5,
                  beta = rep(0,p))
parameter$beta[1:4] <- c(1,2,-1,-2)
#=======================================#
t <- seq(60,120, length.out = 10) #time values
set.seed(123)

G <- 40 ; ng = 4
link = m

dt <- create_JLS_HD_data(G, ng, t, m, link, parameter)

var.true <- dt$var.true
a <- var.true$a ; var.true$a <- NULL #a fixé (et retiré des variables latentes)
S.data <- dt$survival
U_full <- dt$U

Y <- do.call(get_obs, var.true) + rnorm(n, 0, sqrt(parameter$sigma2))
S.data.time <- S.data$obs
S.data.time.log.sum <- sum(log(S.data.time))
model <- SAEM_model( 
  function(sigma2, ...) -n/(2*sigma2),
  function(phi1, phi2, phi3, ...) mean((Y - get_obs(phi1, phi2, phi3) )^2 ), 'sigma2',
  
  # === Variable Latente === #
  latent_vars = list(
    # === Non linear model === #
    latent_variable('phi', dim = G, size = 3, prior = list(mean = 'mu', variance = 'omega2'),
                    add_on = c('zeta(phi1 = phi1, phi2 = phi2, phi3 = phi3, ...)' )),
    
    # === S.data model === #
    latent_variable('b', prior = list(mean = 'barb', variance.hyper = 'sigma2_b'),
                    add_on = c('zeta(b = b, ...) +',
                               'sum(h$eval(b = b, ..., i = c(1,2)))' )),
    latent_variable('alpha', prior = list(mean = 'baralpha', variance.hyper = 'sigma2_alpha'),
                    add_on = c('zeta(alpha = alpha, ...) +',
                               'alpha*h$eval(alpha = alpha,..., i = 3)'))
  ),

  # === Paramètre de regression === #
  regression.parameter = list(
    regression_parameter('beta', 1, function(...) SPGD(5, theta0 = beta,
                                                      step = 0.05, lambda = 1/sqrt(N),
                                                      normalized.grad = T,
                                                      zeta.der.B, N, zeta.B, 
                                                      Z$alpha,  Z$phi1, Z$phi2, Z$phi3,Z$b) )
  )
)
# ---  Initialisation des paramètres --- #
parameter0 <- parameter %>% sapply(function(x) x* runif(1, 1.1,1.4))
#===============================================#
load.SAEM(model)
U <-  U_full
oracle <- maximisation(1, do.call(S$eval, var.true), parameter, var.true)
#==============================================================================#

init.options <- list(x0 = list(phi = c(1,80,4), b = parameter0$barb, alpha = parameter0$baralpha), 
                     sd = list(phi = c(.05, 1.5, .5), b = 1, alpha = .1) )

SAEM.options <- list(niter = 300, sim.iter = 5, burnin = 190, 
                adptative.sd = 0.6)

saem

p <- 4
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 32min 38sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.0000 0.5000
Estimated value 0.0025 0.9034 89.9207 5.0293 0.0039 35.3106 0.6393 28.7713 0.8068
Rrmse 0.0008 0.0002 0.0004 0.0043 0.0194 0.0169 0.0798 0.0410 0.6135 
p <- 20
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 36min 24sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.0000 0.5000
Estimated value 0.0025 0.9045 89.9446 5.0482 0.0039 35.5880 0.6418 33.5285 0.6676
Rrmse 0.0023 0.0015 0.0001 0.0081 0.0152 0.0092 0.0762 0.1176 0.3352 
p <- 100
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 42min 04sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.0000 0.5000
Estimated value 0.0025 0.9034 89.9501 5.0219 0.0039 35.0544 0.6071 36.9021 2.2173
Rrmse 0.0012 0.0002 0.0001 0.0028 0.0199 0.0240 0.1262 0.2301 3.4346 
p <- 500
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 37min 12sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.0000 0.5000
Estimated value 0.0025 0.9028 89.9222 5.0191 0.0038 34.9631 0.6428 25.2145 1.1898
Rrmse 0.0022 0.0004 0.0004 0.0023 0.0431 0.0266 0.0748 0.1595 1.3797 
p <- 1000
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
saveRDS(res, paste0(params$rds_filename, '_', p, '.rds'))
## [1] "SAEM execution time = 00h 34min 15sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.0000 0.5000
Estimated value 0.0025 0.9030 89.9085 5.0218 0.0039 35.1083 0.6076 22.1310 0.5723
Rrmse 0.0006 0.0003 0.0005 0.0028 0.0100 0.0225 0.1255 0.2623 0.1447